Histogram Monte Carlo Simulation of the 
Geometrically Frustrated XY 
Antiferromagnet with Biquadratic Exchange 

M.Zukovic^*, T.Idogaki^ and K.Takeda^-^ 

^ Department of Applied Quantum Physics, Kyushu University 
^ Institute of Environmental Systems, Kyushu University 

Abstract. Histogram Monte Carlo simulation is used to investigate effects of biquadratic 
exchange J2 on phase transitions of a 3D classical XY antiferromagnet with frustration in- 
duced by the antiferromagnetic exchange Ji and the stacked triangular lattice geometry. 
The biquadratic exchange is considered negative (antiferroquadrupolar) within the triangu- 
lar planes and positive (ferroquadrupolar) between the planes. The phase diagram obtained 
features a variety of interesting phenomena arising from the presence of both the biquadratic 
exchange and frustration. In a strong biquadratic exchange limit (|t/i|/|t/2| < 0.25), the an- 
tiferroquadrupolar phase transition which is of second order is followed by the antiferromag- 
netic one which can be either first or second order. The separate antiferroquadrupolar and 
antiferromagnetic second-order transitions are found to belong to the chiral XY and Ising 
universality classes, respectively. If the biquadratic exchange is reduced both transitions are 
found to be first order and occur simultaneously in a wide region of |Ji|/| J2|- However, if 
IJ2I — > the transition changes to the second-order one with the chiral universality class 
critical behavior. 
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I. Introduction 



The problem of biquadratic or generally higher-order exchange interactions in systems 
with Heisenberg symmetry has been addressed in several mean field approximation (MFA) 
studiesi-^, by high-temperature series expansion (HTSE) calculations^, as well as within a 
framework of some other approximative schemes^i^. It has been shown that such interactions 
can induce various interesting properties such as tricritical and triple points, quadrupole or- 
dering, separate dipole and quadrupole phase transitions etc. Much less attention, however, 
has been paid to this problem on systems with XY spin symmetry. Chen et.al—'^ calcu- 
lated transition temperatures and the susceptibility critical indices for an XY ferromagnet 
with biquadratic exchange on cubic lattices by the HTSE method for limited region of Ji/ J2- 
However, the rigorous proof of the existence of dipole long-range order (DLRO), correspond- 
ing to the ferromagnetic directional arrangement of spins, and quadrupole long-range order 
(QLRO), representing an axially ordered state in which spins can point in either direction 
along the axis of ordering, at finite temperature on the classical bilinear-biquadratic exchange 
model has only recently been provided independently by Tanaka and Idogaki^, and Campbell 
and Chayes^. Very recently we have considered the XY model with the bilinear-biquadratic 
exchange Hamiltonian on a simple cubic^^ and hexagonal (stacked triangular) lattices, and 
performed a finite-size scaling (FSS) analysis in order to investigate critical properties of the 
considered systems via Standard Monte Carlo (SMC) and Histogram Monte Carlo (HMC) 
simulations. 

So far, however, to our best knowledge there has been no investigation of the effect of 
the biquadratic exchange on an XY model with frustrated/competing exchange interaction. 
In this paper we present systematic investigations of the role of the biquadratic exchange in 
phase transitions of the geometrically frustrated XY antiferromagnet on stacked triangular 
lattice (STL). This model has been argued to possess some unique properties such as novel 
chiral universality class critical behavior— li^, but many more remarkable features have been 
observed when the effects of external magnetic field^^ and next-nearest neighbors^ were 
considered. In the present work, the effect of the biquadratic exchange is also found to bring 
about variety of interesting phenomena, such as regions of first order transitions, separate 
magnetic and quadrupolar ordering, transitions of different universality classes, etc. 
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II. Model and computation details 



We consider the XY model, described by the Hamiltonian 



H — —Ji Si ■ Sj — J2 ^^{Si ■ Sk)^ — J2 ^(-^i ■ "^0^ ' 



(1) 



{i,k) {i,l) 



where Si = {Sf,Sf) is a two-dimensional unit vector at the ith lattice site and the sums 
{i, k) and {i, I) run over all nearest neighbors (NN), NN in the xy-p\ane, and NN in the 
stacking 2;-axis direction, respectively. We consider the bilinear exchange interaction Ji < 0, 
the biquadratic intra- plane and inter-plane exchange interactions < and > 0, re- 
spectively, with |J^| = I j|| = IJ2I. 

Assuming periodic boundary condition, spin systems of the linear lattice sizes L — 12, 
18, 24 and 30 are first used in SMC simulations. For a fixed value of the exchange ratio 
|<>^i|/|^2|) wc start the simulation process at low (high) temperatures from an antiferromag- 
netic/random (random) initial configuration and gradually raise (lower) temperature. These 
heating-cooling loops serve to check possible hysteresis, accompanying first-order transitions. 
As we move in (| Ji|/| J2I, ^b^/|<-^2|) space, we use the last spin configuration as an input for 
calculation at the next point. We sweep through the spins in sequence and updating follows 
a Metropolis dynamics. In the updating process, the new direction of spin in the spin flip 
is selected completely at random, without any limitations by a maximum angle of spin ro- 
tation or allowed discrete set of resulting angle values. Thermal averages at this stage are 
calculated using at most 1 x 10^ Monte Carlo steps per spin (MCS/s) after thermalizing 
over another 0.5 x 10^ MCS/s. We calculate the system internal energy E and some other 
physical quantities defined as follows: the specific heat per site c 



where is the ath sublatticc-magnctization vector (note that the present model has six 
equivalent magnetic sublattices) , given by 



(2) 



the dipole LRO (DLRO) parameter m. 




(3) 




(4) 
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the quadrupole LRO (QLRO) parameter g, 



N N 



where 



\ i i 

the chiral LRO (CHLRO) parameter 



(5) 



(6) 



K, — 



N N 



(7) 



where the summation runs over all upward triangles on the triangular layer and Kp represents 
a local chirality at each elementary triangular plaquette, defined by 



3^3 



Y^[Si X Sj]^ = ^^[sin((/j2 - (Pi) + sin((^3 - (^2) + sin((^i - (^3)] , (8) 



(hj) 



where the summation runs over the three directed bonds surrounding each plaquette, p, and 
(pi represents the ith spin angle. Kp is an Ising-like quantity representing the sign of rotation 
of the spins along the three sides of each plaquette. Further, the following quantities which 
are functions of the parameter O (= M, Q, K) are defined: the susceptibility per site xo 

im {cm 



Xo 



NkeT 



(9) 



the logarithmic derivatives of (O) and (O^) with respect to P — l/ksT 
the fourth-order long-range order cumulant U (Binder parameter) 



Upl- 
and the fourth-order energy cumulant V 

V = 1 



3(02)2 ' 



(10) 

(11) 

(12) 



3{E 



2\2 



(13) 
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The above quantities are useful for localization of a transition as well as for determination of 
its nature. For example, first-order transitions usually manifest themselves by discontinuities 
in the order parameter and energy, and hysteresis when cooling and heating. If transition is 
second order, it can be localized approximately by the xo peak position or more precisely 
by the intersection of the fourth-order LRO (or energy) cumulants curves for different L. 

In order to increase precision and reliability of the obtained information, as well as to re- 
trieve some additional information which could not be extracted from the SMC calculations, 
we further perform HMC calculations, developed by Ferrenberg and Swendseni2.iil^ at the 
estimated transition temperatures for each lattice size. Here, 2 x 10^ MCS/s are used for 
calculating averages after discarding another 1 x 10^ MCS/s for thermalization. We calculate 
the energy histogram P{E), the order parameters histograms P{0) (O = M, Q, K), as 
well as the physical quantities (121)- (IT^ . Using data from the histograms, one can calculate 
physical quantities at neighboring temperatures, and thus determine the values of extrema 
of various quantities and their locations with high precision for each lattice size. In such a 
way we can obtain quality data for FSS analysis which determines the order of the transition 
and, in the case of a second-order transition, it also allows us to extract critical indices. For 
example, the energy cumulant V exhibits a minimum near critical temperature T^, which 
achieves the value V^* = | in the limit L — oo for a second-order transition, while V* < ^ 
is expected for a first-order transitioni^ii^. Temperature-dependences of a variety of ther- 
modynamic quantities display extrema at the L-dependent transition temperatures, which 
at a second-order transition are known to scale with a lattice size as, for example: 

XO,ma.(i:) OcL^°/^° , (14) 

DlO,macc{L) CX L^/^^ , (15) 
D20,maAL) <X L'/"^ , (16) 

where up and 70 represent the correlation length and susceptibility critical indices, respec- 
tively. In the case of a first-order transition (except for the order parameters), they display 
a volume-dependent scaling, oc L^. The simulations were performed on the vector super- 
computer FUJITSU VPP700/56. 

III. Chirality on frustrated quadrupoles 
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It has been known for some time that the frustrated spin system on triangular lattice 
possesses the chirality n as defined in Eqns. (17|^ ^^. Due to the chirality the system has two- 
fold degeneracy of the ground state (k = +1 and k = —1), resulting in the structure with 
spins arranged on plaquettes with turn angles +120° and —120°, respectively (Fig.l(a)). 
A minimum energy condition is realized by an arrangement in which the + and — pla- 
quettes alternate, producing long-range chiral order at low temperatures. Such a system 
has been argued to belong to a nonstandard universality class linked to the two-fold chiral 
degeneracy inherent to the 120° ordered spin structure^ii^, the critical behavior of which 
is characterized by critical indices, different from those for non-frustrated systems with the 
same spin symmetry. Since the present Hamiltonian includes both bilinear and biquadratic 
terms, let us take a closer look at the opposite side of the exchange ratio spectrum and 
investigate critical behavior of the system with only biquadratic exchange interaction, i.e. 
the case of Ji = 0. If < (the sign of j| is irrelevant in the present consideration) 
the quadrupolar system is frustrated due to the triangular lattice geometry, resulting in a 
non-collinear ground state. The non-collinear ground state arrangement resembles the 120° 
structure of the antiferromagnetic system, however, here, the spins can point in either di- 
rection within the given axis (for illustration see the snapshots in Fig. 12). As far as the 
chirality n is concerned, such a system has four-fold degeneracy in the ground state of each 
plaquette (k^ = ±1,±|), resulting in the structure with four possible turn angles between 
two neighboring spins ±120°, ±60°. However, there is no energetically favorable arrange- 
ment among the four kinds of plaquettes and, hence, the plaquettes do not order even at low 
temperatures. Nevertheless, even for such a system we can define the quantity analogous 
to the chirality of the antiferromagnetic system (let us call it the quadrupolar chirality) if 
we consider instead of spins their axes and turn angles between the axes, which are again 
±120°. If we define the local quadrupolar chirality as 



concerning such defined quadrupolar chirality, the system will have two-fold degeneracy of 
the ground state {k'^ = —1 and = ±1, corresponding to turn angles ±120° and —120°, 
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= ^--/=[sin2(v?2 - V^i) ± sin2(v93 - (^2) + sin2(v9i - (ps)] , 
and the quadrupolar chirality LRO parameter (QCHLRO) n'^ as 



(17) 




(18) 



respectively (Fig.l(b))), and the situation will much resemble the one for the antiferromag- 
netic system with the chirality k. Furthermore, in analogy with the chirality k which is 
believed to order along with spins, here, the quadrupolar chirality nfl is expected to show 
LRO simultaneously with quadrupoles. 

IV. FSS analysis and phase diagram 

We first consider the case of J2 = 0. To determine the order of the transition we ana- 
lyze the scaling behavior of the minimal value of the energy cumulant V at the transition 
temperature. As shown in Fig.2, V tends to the value of 2/3, as expected for a second- 
order transition, and the slope 2.39 means that V is not volume dependent. Also, observing 
the energy and LRO parameters distribution histograms (not shown), no bimodal distribu- 
tion, which would signal a first-order transition, is found. Hence, both spin and chirality 
ordering transitions seem to be clearly of second order. The transition temperature, cal- 
culated from the intersection of the Binder parameter curves for different L, is estimated 
to ksTc/lJil = 1.4580 ± 0.0005, in agreement with the values quoted in Refs.— li^. The 
chirality transition temperature ksT^/lJil = 1.4590 ± 0.0013, similarly as in Ref.— , seems 
to be slightly higher than the spin ordering temperature but the two values cannot be dis- 
tinguished beyond the error bar and, hence, we assume they are the same. The spin and 
chirality critical indices calculated from the scaling relations ( IT4|) - (fT6|) take the following 
values: um = 0.52 ± 0.03, 7m = 1-08 ± 0.08 and uk = 0.55 ± 0.01, -fx = 0.81 ± 0.03^^, 
respectively (Figs. 3,4). Also the values of the critical indices are in fair agreement with the 
two previous studies^^ii^, however, as far as the universality class is concerned the situation 
here is not so straightforward and will be discussed later. 

The order of the transitions changes, however, when even a comparatively weak bi- 
quadratic exchange interaction is introduced. Although it is very hard to observe the typical 
first-order behavior for small values of | J2I, if the lattice sizes are taken sufficiently large the 
signs of the discontinuous transition show up. This is seen in Fig.5 in which the bimodal 
(double-peak) energy distribution becomes clearly recognizable if L > 30, for the case of 
|^2|/|</i| = |- As I J2 1 is increased, the first-order features of the transition are becoming 
more and more apparent. Fig. 6 shows clearly bimodal energy distribution histograms for 
|<^i|/|<^2| = 1-3, in which the dip between the peaks is observable already at smaller L, quite 
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rapidly approaching zero as L is increased, indicating discontinuous behaviour of the energy 
at a rather strong first-order transition. Although we do not show it here, similar double 
peaks can also be observed in the histograms of each LRO parameter. 

The transition remains first order and simultaneous for dipole, quadrupole and chiralities 
ordering until fairly small values of | Ji|/| J2|- Below | Ji|/| J2I — 0.25, however, quadrupoles 
order separately at temperatures higher than those for dipole ordering. Thus the phase 
boundary branches and a new middle phase of axial quadrupole long-range order (QLRO) 
without magnetic dipole ordering opens between the paramagnetic and DLRO phases. This 
phase broadens as | Ji | / 1 J2 1 decreases, since the QLRO branch is little sensitive to the | Ji | / 1 J2 1 
ratio variation and levels off, while the DLRO branch turns down approaching the point 
(|-^i|/|-'2|, ^^2^/1-^21) = (0,0). This means that the ground state is always magnetic as long 
as there is a finite dipole exchange interaction. In Fig. 7 we present the temperature variation 
of the DLRO, QLRO, CHLRO and QCHLRO parameters m, q, k and k*, respectively, at 
|<^i|/|<^2| = 0.15. We can see that quadrupoles order before dipoles, forming a fairly broad 
region of QLRO without DLRO. On the other hand, the chirality and quadrupole chirality 
seem to order simultaneously with dipoles and quadrupoles, respectively. The QLRO transi- 
tion is apparently second order down to | Ji|/| J2I = and the critical indices take the values 
VQ = 0.50 ± 0.03, 7Q = 1.09 ± 0.08 at | Ji|/| J2I = 0.15 (Fig.8) and vq = 0.520 ± 0.003 and 
7q = 1.072 ± 0.009 at Ji = 0. In the case of Ji = 0, the QLRO transition temperature is 
located as fc^Tg/l J2I = 0.729±0.002. On the other hand, in the case of the DLRO transition, 
the first order seems to persist even after the QLRO and DLRO boundaries separate for a 
small range of the exchange ratio values just below the splitting point. This is clearly seen in 
Fig. 9 from the distribution diagrams of the DLRO and QLRO parameters. Although at first 
glance it seems that both transitions occur at the same temperature and are of first order, a 
closer look reveals that while the bimodal distribution of the DLRO parameter is between the 
disordered and ordered states, the bimodal distribution of the QLRO parameter is between 
two ordered states of different finite QLRO parameter values. Therefore, here, the QLRO 
parameter only shows a discontinuity within the QLRO region, rather than paramagnetic- 
QLRO transition. The first-order DLRO transition changes to the second-order one upon 
further lowering of | Ji|/| J2|. This is seen from the finite-size scahng analysis of the HMC 
data for | Ji|/| J2I = 0.15 (Fig. 10). The slopes apparently indicate the second-order character 
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of the transition with the critical indices um = 0.63 ± 0.02, 7m = 1.25 ± 0.04. The resulting 
phase diagram is drawn in Fig. 11 and some relevant numerical results listed in Table 1. 

V. Summary and discussion 

We studied effects of the biquadratic exchange on the phase diagram of the frustrated 
classical XY antiferromagnet on STL. This study, which to our best knowledge is first for 
the studied system, covered most of the significant phenomena induced by the presence of 
the biquadratic exchange, and present a fairly compact picture of the role of this higher- 
order exchange interaction on the critical behaviour of the system considered. We obtained 
the phase diagram with two ordered phases: in the region where the bilinear exchange is 
dominant there is a single phase transition to the DLRO phase, which is second order at 
J2 = 0, but changes to a first-order one upon adding of a rather small amount of biquadratic 
exchange. In the region of small |t/i|/|J2| the phase boundary splits into the QLRO transi- 
tion line at higher temperatures and the DLRO transition line at lower temperatures, which 
are second order, and partly first and partly second order, respectively. 

From our qualitative and quantitative evaluations we found out that not only the order 
of the transitions in different regions of the IJ1I/IJ2I parameter is not the same, but also 
the sets of the critical indices obtained in different regions of the second-order transition 
are different for seemingly the same kind of transition while almost identical for different 
kinds of transition. Let us first discuss the problem of the order of the transition. The 
second-order transition at J2 = is in agreement with the previous MC studiesi^ii^ but 
in contradiction with the renormalization group study22^ which predicts a clear first-order 
transition. At finite J2, the first-order transition observed in the region of the paramagnetic- 
DLRO transition has also been observed in the case of a ferromagnet with Ji > 0, > 
and J2 > 0, however, only in a quite narrow region of J1/J2 ^ (0.33,0.55)^^. We believe 
that the mechanism responsible for this transition in the present case is similar to that in 
the case of the ferromagnet, i.e., it could result from a kind of tension between the bilin- 
ear and biquadratic exchange interactions, which in the present case seems to be enhanced 
by the presence of the frustration and consequently causing broadening of the first-order 
transition region. Namely, while the decreasing bilinear exchange drives the transition tem- 
perature down to the lower values, the biquadratic exchange does not follow this tendency 
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and rather prevents the ordering temperature from rapid decrease. This tendency is clearly 
seen from the phase diagram both in the region of separate transitions, where Tg does not 
vary much with decreasing |Ji|/| J2I, as well as in the region of simultaneous ordering, where 
the transition temperature is apparently enhanced by the presence of the biquadratic ex- 
change (the case of absent biquadratic exchange is represented by the dash-dot straight hne 
in (|Ji| — fcfiTc) parameter space). Put differently, quadrupoles would prefer ordering at 
higher temperatures but as long as there is a single transition they are prevented to do so by 
too low bilinear exchange, and order occurs only if the temperature is lowered still further. 
This "frustration" results in a first-order transition when the strength of the quadrupole 
ordering prevails and frustrated quadrupoles order abruptly along with dipoles. However, 
when I J2 1 reaches high values the frustration becomes too high for the two kinds of ordering 
to occur simultaneously and they separate. In order to understand the first-order DLRO 
transition and QLRO parameter discontinuity in the region just below the point of the sep- 
aration, we analyzed snapshots (not shown) for | Ji|/| J2I = 0.25 just before the DLRO sets 
in. In the snapshots we could observe fairly large clusters of antiferromagnetically ordered 
spins along the stacking direction, which is non-frustrated and in which spins seem to or- 
der more easily than within frustrated planes (Note that in the case of the non-frustrated 
parallel (ferromagnetic) ordering the transition temperature is roughly twice as higher as 
in the present case^^). These clusters reorient at the transition as a whole, and such a way 
may produce discontinuities in the order parameter and internal energy i.e., a first-order 
transition. Besides those clusters, we could also observe smaller intra-plane clusters of spins 
the axes of which show local parallel ordering. At the DLRO transition, the spins in these 
clusters (and also their axes) reorient into the 120° spin structure, which may result in the 
small discontinuity of the QLRO parameter, seen in Fig.9. The separate QLRO is appar- 
ently second order, in agreement with the mapping arguments of Carmesin^i. 

Now, let us address the problem of the critical indices in the case of a second-order tran- 
sition. For the case of J2 = 0, there has been an argument about the universality class of 
the critical behavior of such a system. Kawamura claimed that it should display a nonstan- 
dard (chiral) universality class behavior due to a two-fold discrete degeneracy Z2 associated 
with the two chiral states, with the novel indices {om = 0.34 ± 0.06, Pm = 0.253 ± 0.01, 
7m = 1.13±0.05 and um = 0.54 ± 0.02)^, while Plumer et alM maintained that there is no 
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new universality class and that the indices take the mean-field tricritical values (om = |, 
I^M = \, 1m = 1 and pm = \)- As we can see, the DLRO critical indices obtained from 
our calculations vm = 0.52 ± 0.03, 7m = 1.08 ± 0.08 are somewhere between those from 
Refs.— >^ and, considering the error estimates, could be interpreted for support of either of 
the theories. Although we can make no definite conclusion based on the values of the indices 
themselves, we believe that the former interpretation is more favorable. Indeed, looking at 
the critical indices of the separate QLRO transitions we can see that they are strikingly 
similar to those for the case of J2 = (and seem to be such along the whole paramagnetic- 
QLRO boundary). These indices can hardly be interpreted as the mean-field tricritical ones 
and the theory of the same universality class critical behavior of quadrupoles (Ji = 0) 
and dipoles (J2 = 0), based on mapping and quantitative analysis^^, would rather strongly 
suggest that both cases show the chiral universality class behavior. As far as the separate 
DLRO transition is concerned, in the second-order transition region we obtained the critical 
indices um = 0.63 ± 0.02, 'Jm = 1-25 ± 0.04, which are quite different from those for the 
DLRO transition at J2 = 0. The reason is that in this case the transition has an Ising-like 
character and, hence, the indices take on the Ising universality class values (i/^''*"^ = 0.629, 
^ising ^ 1,23924). The situation is illustrated in Figs.l2(a,b). In the QLRO (and no DLRO) 
region, there is an axial quadrupole ordering on each of the three sublattices (Fig. 12(a)) 
and only upon further lowering of the temperature the system reaches the QLRO-I-DLRO 
phase in which the Ising-like directional dipole ordering within the given axis in each sublat- 
tice takes place (Fig. 12(b)). Therefore, here, the only difference from the Ising case is that 
dipoles can order along any of the three axes, not only the z-axis. 

Our further intention is to perform similar simulations on the STL antiferromagnet for 
some other interesting cases, like: J2 < 0, j| < 0; > 0, j| > 0; > 0, j| < 0. 
Besides the geometrical frustration, such spin systems will possess additional frustration 
arising from the bilinear and biquadratic exchanges competing in the stacking direction, 
intra-plane direction, and both stacking and intra-plane directions, respectively. 
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(b) 

FIG. 1: Two degenerate ground states, +120° and —120° structures on (a) spin and (b) quadrupole 
plaquettes. Signs + and — denote the sign of (a) chirality and (b) quadrupole chirality of the 
elementary triangles. Spins and quadrupoles are numbered counterclockwise, corresponding to the 
definitions ([8]) and pT)) . 



0.003 



0.002 



0.001 



^ ' ' ' ' ' 

0.0005 0.001 0.0015 0.002 0.0025 

L-2.39 

FIG. 2: Scaling of the energy cumulant minima at J2 = 0. The values extrapolated to L ^ 00 
approach the value F* = | and do not scale with volume, as it should be in the case of a second- 
order transition. 



TABLE I: Critical indices and transition temperatures for quadrupole, dipole, and chiral ordering, 
respectively. 
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FIG. 3: Scaling behaviour of the maxima of the susceptibility XM,max corresponding to the pa- 
rameter M and logarithmic derivatives of its first and second moments Dim, max and D2M,maxj 
respectively, in In-ln plot, for J2 = 0. The slopes yield values of I/i^m for DiM,max, D2M,max and 
Im/^m for 
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FIG. 4: The same dependence as in Fig.3, with the parameter K considered instead of M. 
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E/N 



FIG. 5: Energy distribution at the size-dependent transition temperatures Tc{L) for various lattice 
sizes and | J2|/|'^i| = ^- The bimodal distribution signaUng a first-order transition can only be seen 
at L > 30. 
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FIG. 6: Energy distribution at Tc{L) for |Ji|/|J2| = 1-3. Doublc-pcaked structure with deepening 
barrier between the two energy states with increasing lattice size indicates a first-order transition. 
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FIG. 7: Temperature variation of the DLRO, QLRO, CHLRO and QCHLRO parameters m, q, n 
and K*, respectively, for |Ji|/|J2| = 0.15 and L = 12. 
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FIG. 8: Scaling behaviour of the maxima of the susceptibility XQ,max and logarithmic derivatives 
of the parameter Q and its second moment -DiQ,maa; and D2Q^max, respectively, in In-ln plot, for 
|'^i|/|'/2| = 0.15. The slopes yield values of 1/z/q for Dig^max, D2Q,max and jq/uq for XQ,max- 
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0.025 




M/N, Q/N 

FIG. 9: Distribution histograms P[M) and P{Q) of DLRO and QLRO parameters, respectively, at 
Tc{L) for IJ1I/IJ2I = 0.25. The bimodal distributions of the DLRO and QLRO parameters signal 
a first-order disorder-DLRO transition and a jump between two finite values of QLRO parameter, 
respectively (see text). 



20 



s 



Q 



^■■•ln(DiM,max)' slope 


= 1.59 


<>. . . ln(D,„ „„J, slope 


= 1.58 


0...1n(XM,max). slope = 


1.98 







2.4 



2.8 



3.2 



3.6 



ln(L) 



FIG. 10: Scaling behaviour of the maxima of the susceptibiHty XM,max and logarithmic derivatives 
of the DLRO parameter and its second moment Dim, max and D2M,maxj respectively, in In-ln plot, 
for I J2I = 0.15. The slopes yield values of 1/i/m for DiM,max, D2M,max and ^m/^^m for XM,max- 
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FIG. 11: Phase diagram in (| Ji|/| J2I, ^B^c/I J2I) space. The paramagnetic (P), antiferroquadrupo- 
lar (AFQ), and antiferromagnetic (AFM) regions correspond to the phases in which both dipoles 

and quadrupolcs arc disordered, only quadrupolcs arc ordered, and both dipoles and r^uadrupoles 
are ordered, respectively. The solid and dashed lines correspond to second- and first-order transi- 
tions, respectively. 
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FIG. 12: Spin configuration snapshots of the system for |Ji|/|J2| = 0.05 in (a) QLRO phase 
(fcsT/lJal = 0.3) and (b) DLRO phase (fcsT/lJal = 0.001). 



23 



